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ABSTRACT 

Simulations show eccentric disks (m = 1 modes) forming around quasi-Keplerian potentials, 
a topic of interest for fueling quasars, forming super-massive BHs, planet formation and mi- 
gration, explaining the origin and properties of nuclear eccentric stellar disks like that in M31, 
and driving the formation of the obscuring AGN torus. We consider the global, linear normal 
m = 1 modes in collisionless disks, without the restriction that the disk mass be negligible 
relative to the central (Keplerian) mass. We derive their structure and key resonance features, 
and show how they arise, propagate inwards, and drive both inflow/outflow and eccentricities 
in the disk. We compare with hydrodynamic simulations of such disks around a super-massive 
BH, with star formation, gas cooling, and feedback. We derive the dependence of the normal 
mode structure on disk structure, mass profiles, and thickness, and mode pattern speeds and 
growth rates. We show that, if the disk at some radii has mass of > 10% the central point 
mass, the modes are linearly unstable and are self-generating. They arise as "fast modes" with 
pattern speed of order the local angular velocity at these radii. The characteristic global nor- 
mal modes have pattern speeds comparable to the linear growth rate, of order (GMqRq 3 ) 1 ^ 2 , 
where Mq is the central mass and Ro is the radius where the enclosed disk mass ~ Mq. They 
propagate inwards by exciting eccentricities towards smaller and smaller radii, until at small 
radii these are "slow modes." With moderate amplitude, the global normal modes can lead to 
shocks and significant gas inflows at near-Eddington rates at all radii inside several ~ R®. 

Key words: quasars: general — galaxies: active — galaxies: evolution — galaxies: nuclei — 
cosmology: theory 



1 INTRODUCTION 

Perturbations to normally circular orbits in nearly-Keplerian po- 
tentials are of fundamental interest for a variety of topics in as- 
trophysics. For example, questions related to the fueling of super- 
massive black holes (BHs), their accretion disks, and the dynam- 
ics of nearby systems, the formation and fueling of protostars, and 
the behavior of protoplanetary disks and planets around stars or 
rings and moons around planets. Particularly interesting are pertur- 
bations with azimuthal wavenumber m = 1 (amplitude oc cos<j!>), 
which can manifest as eccentric orbits or disks, lopsided or slosh- 
ing modes, or one-armed spirals. It is easy to see why: the response 
of a nearly circular orbit to a weak perturbation, to leading or- 
der, scales with 1 /[k 2 — m(Q — Q p ) ], where Q is the orbital fre- 
quency, k is the epicyclic frequency, and Q p is the characteristic 
frequency (precession rate) of the perturbation. In a Keplerian po- 
tential, k — Q oc r _3//2 , so (since O p is finite) for any continuous 
system this scales at small radii as ~ 1/(1 — m) Q, 2 . For general m, 
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this vanishes, but for m = 1 the leading terms cancel and there is 
a strong resonant response. Physically, this 1:1 resonance between 
radial and azimuthal frequencies is related to the fact that ellipti- 
cal orbits in a Keplerian potential are closed and do not precess. 
As a consequence, the eccentricity distribution and mode behavior 
in such a disk can be determined by collective effects in the disk, 
even where these collective effects are very weak compared to the 
gravity of the central object. 

Recently, for example, |Hopkins & Quataert (2009) have 
shown that the formation of lopsided, eccentric disks within the BH 
radius of influence is a ubiquitous feature in hydrodynamic simu- 
lations of massive gas inflows in galaxies, and that such disks can 
efficiently drive gas angular momentum loss and power BH accre- 
tion rates of up to ~ IOMq yr _ 1 . The co-existence of gas and stars 
is critical for the large inflow rates seen; the torques on the gas are 
dominated by the mode in the collisionless portion of the disk. And 
the stellar relics of these disks bear a remarkable similarity to nu- 
clear disks observed on < lOpc scales around nearby supermassive 
BHs, particularly the well-studied case at the center of M31 jLauer] 
|et al.| 19 93 1, whose origin has been mysterious. The inflow and out- 
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flow regulated by these modes in such simulations also determines 
the nature of the galaxy mass profile on scales within < 10— 100 pc. 
As such, it is also particularly interesting to understand whether or 
not such modes could arise generically. There are many candidate 
nuclear disks in such systems JLauer et al.| 1996 2005] [Houghton| 
et al.|2006||Thatte et al. 2000; Debattista et al.|2006||Afanasiev & 



Sil'chenko 2002, Seth et al. 2010). 



There is considerable literature discussing the mode structure, 
pattern speeds, and evolution of general self-gravitating disk in- 
stabilities (see e.g. Lin et al. 1969; Goldreich & Tremaine 1978 
1979; Toomre 1969 1977). But these m = 1 modes are less well- 
understood, especially in collisionless (stellar or planetary) disks. 
For example, there remains considerable debate regarding the sta- 
bility of such modes (e. g.|Tremaine|2001 [|SaTow & Statler|2001 [[jiT] 



|cobs & Selrwood|2001| r rouma|2002) . These works describe many 
interesting behaviors of m = 1 modes in a disk in a nearly Keple- 
rian potential, but their conclusions rely on specific assumptions. 
Tremaine ( 2001 ) show that such modes are linearly stable, but only 
in the limit where Mj <C Mbh at all radii. And much of the insight 
from the study of protostellar and planetary disks focuses on ei- 
ther the same limit (for planetary disks) or the opposite limit (for 
early-stage protostellar disks) in which the system is really just a 
self-gravitating disk (see |Laughlin & B odenheimer 1994; Laugh- 
|lin & R ozyczka 1996, and references therein) But in simulations 
or observations of galactic nuclei JLevine et al.|2008 ; Escala|2007| 



Mayer et al.|2007||Hopkins et al. 2010; Hopkins & Quataert 2009) 



and protostellar evolution 



Laug hlin & Bod enheimer 1994 , [Bate] 



|et al.| [l995; Nel son et aT]|l998| >, much of the interesting behav- 
ior involves disks over a range of radii where the disk is not com- 
pletely negligible in mass. The mode structure outlined in Tremaine 
(2001 ) is used in Papaloizou (2002) to estimate the effects on a one 
or two low-mass planet system, but assuming mode stability and 
under similar mass and radius restrictions (and not allowing for a 
large collisionless disk component). [Adams et aL|( |1989| > and |Os-| 
triker et al. ( 1992); Shu et alT|( |1990| reach opposite conclusions to 
Tremaine (2001 1, but they focus on "fast modes" (Q p ~ Q.) in the 
regime Mj ~ Mbh (or Mj ~ M,„ ). and consider only pure fluid disks 
with a hard and reflecting outer edge, in which case the resulting 
mode growth rates depend sensitively on the structure and nature 
of reflection at the disk edge (and can be dramatically modified - 
in fact, completely eliminated - allowing any flow "through" or 
outside the edge, the most probable configuration). In general, col- 
lisionless disks can have very different mode structure from fluid 
systems. Christodoulou & Narayan ( 1992) consider instabilities of 
narrow torus annuli (as opposed to disks) but with similar restric- 
tions. Also, because the mode growth is sensitive to motion of the 
center of mass, most numerical studies to date have lacked the res- 
olution to determine whether mode growth is real or artificial in the 
limit where the disk mass is smaller than the BH/star mass (see the 
discussion in Nelson et al. 1998). 

Moreover, these studies have largely been restricted to very 
specific mass profiles (e.g. the Kuz'min disk, with £ — ^constant at 
small radii) and/or systems with sharp "edges," where in fact the 
profiles seen in interesting phases in simulations and in e.g. star- 
forming systems and the centers of real "cusp" elliptical galaxies 
resemble a range of power-law slopes with smooth declines at large 
radii and significant variation in profile shape ( Gebhar dt et al.| 1996| 
|Hopkins et al.|2009| [Hopkins & Hernquist|2010) . ~ 

As such, the origin of these modes, where observed, and many 
of their properties have remained ambiguous. Given claims of the 
stability of such modes, alternative suggestions for their origin have 
ranged from their being induced by an external collision/passage 



(in galactic nuclei, from e.g. a nuclear star cluster; see Sambhus 
& Sridhar 2002 1, to their being excited by substantial populations 
of stars on retrograde orbits (Touma 2002 ). But these explanations 
pertain only to very specific systems or regimes, and do not ex- 
plain the ubiquity of such modes in astrophysical systems. Showing 
that they can in fact be self-generating via gravitational instability 
would have profound implications. 

Moreover, in systems with star or planet formation as an ongo- 
ing process, these modes will evolve in time, as quantities such as 
the disk gas fraction, mass profile, dispersion profile, and total disk 
mass are affected by these processes. Similarly, the modes in the 
collisionless disk can drive very strong inflows and outflows in gas, 
changing the properties of t he host disks in turn (see | Nelson et al.| 



[T9^[Bournaud et al.|2005||Hopkins & Quataert|2009| >. Therefore 
it is of great importance to both survey a range of analytic profile 
shapes and dispersion levels, and to compare with simulations that 
can incorporate non-linear effects on the mode evolution. It is also 
particularly interesting to examine the level of inflow generated by 
such modes, if the disks have some gas. 

In this paper we expand upon the investigation of global m=\ 
modes in nearly Keplerian, predominantly collisionless disks. In 
particular, we focus on the question of whether or not such modes 
can, in fact, be unstable and self-generating, and if so how they 
propagate throughout such disks. After defining some terms (§[2), 
we consider modes in the local (WKB) limit (§ |5), which allows 
us to analytically derive approximate stability criteria and discuss 
conditions for efficient mode propagation, in both gaseous and stel- 
lar disks. In §g] we discuss exact numerical solutions which allow 
us to extend this discussion to linear global normal modes, and out- 
line under what conditions these modes are unstable (and what their 
characteristic frequencies are), and how this depends on a variety of 
disk properties including mass profile shape, mass, and disk thick- 
ness. We discuss the structure of such modes, and the conditions 
under which they will drive shocks in collisionless+gaseous sys- 
tems and corresponding inflow/outflow, and the radii over which 
the modes can act. In § [5] we compare to the results from high- 
resolution hydrodynamic simulations which include self-gravity, 
gas cooling, star formation, shocks, inflow/outflow, and non-linear 
effects all not captured in an analytic formulation, and discuss how 
this compares to our analytic insig ht. Finally, in §[6] we summarize 
our conclusions and their implications for astrophysical systems. 



2 DEFINITIONS 

Adopt a cylindrical coordinate frame, (R, <j>, z), and consider an 
initially axisymmetric, thin, planar disk with an arbitrary spherical 
(BH+bulge+halo) component. The coordinate center is located at 
the BH. 

The initial potential in the disk plane can be written $o = 
&o(R), and other properties are defined in standard terms: 



$o = $ (fl) 
<9$ 



v; = r 



_ GM eac (<R) 

dR~ R 

l 



: 'dyn - V c/R 

dR 



9 2 $ 



+ 3S1- 



(2) 
(3) 
(4) 



where V c is the circular velocity, the angular velocity, and k 
the epicyclic frequency. We use c s to denote the sound speed in 
a gaseous disk and a z the vertical dispersion in a stellar disk. 
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Consider a perturbation in the plane and define the perturbed 
surface density field by 



£->£o(fl) + Ei(tf, 



(5) 



We consider a frame rotating with the perturbation pattern speed 
flp. We can represent the perturbed system as a sum of linearly 
independent modes m, Ei = 53mLi ^'»- Since we will focus on 
the behavior of these modes individually, consider (without loss of 
generality), the case of a single mode, 



E,„ = T, a (R) exp{i (m4> — ojt)} 
E n (fi) = Ki?)|E (i?)exp(; f kdR 



(6) 
(7) 



where m is the azimuthal wavenumber, \a\ = \a(R)\ the effective 
mode amplitude, k the radial wavenumber, and the complex uj the 
mode frequency. With these definitions, the mode pattern speed is 
flp = Re(o;)/m, and the linear mode growth rate is 7 = Im(iu) 
(|a| oc exp(+7?)). 

The mode of particular interest is an m — 1 mode, in a quasi- 
Keplerian potential. Simulations in which these inflows form self- 
consistently, including gas inflow, star formation, and feedback, 
typically find that the nuclear stellar and gas distributions are well- 
approximated by power laws on the scales of interest (Hopkins & 
|Quataert|2009^ . We will therefore adopt a true power-law disk as a 
convenient reference model. For such a disk 



E <xR~ 



and it is straightforward to show that 



V c = 2naGER 
Md{< R) =2n(2~rj)~ 1 T,R 2 

r[i-?]r[i±2] 



(8) 

(9) 
(10) 

(11) 



for < rj < 2 (a w r\ for r\ < 1, « 1/(2 — rf) for 1 < r\ < 2). Around 
a central BH this gives 

-(*)+!) 



GM BH 27raGE 



Ro 



2 _ GMm , .27raGEo 

K — r h (i TJ) 



Ro 



-iv+i) 



(12) 



(13) 



We will also discuss finite disks, with a power-law cutoff (mo- 
tivated by analogy to the Kuz'min disk) of the form 



-(3- V )/2 



(14) 



so that E oc R 71 at small radii and E oc R 3 at large radii. The total 
mass in the disk is then 



M d = E a 



(15) 



ay v iy i/2 T[l- V /2] 
¥0) r[(3-7?)/2] ■ 

In this case, Q, (and the potential) must be evaluated numerically, 
but for purposes of interpretation, they can be approximated well 
(exactly at small/large R and with ~ 10% accuracy at R ~ a) by 

-i+ v /2 



^2 ^ GMbh _|_ 27raGE 



Ro 



-(v+i) 



where fi 



-(2-t,)/2 : 



(16) 

(7r 1/2 r[l-r ? /2])/(2Qr[(3-r7)/2]).At small 



radii this is just the normal power-law disk il, at large radii the disk 
portion simply becomes Keplerian, Q, 2 — > G (Mbh +Mj) / r 3 . 

Although we, for convenience, define our canonical model 
with respect to a super-massive BH, many of our conclusions could 
identically be applied to any sufficiently collisionless disk in a 
quasi-Keplerian potential (e.g. protestellar or protoplanetary disks). 
In this case one should simply replace "black hole" with "central 
massive object" (e.g. star), and stellar disk with whatever collision- 
less disk surrounds the system. Because of the scale-free nature of 
the above equations, the rescaling of units is trivial. 



3 THE WKB APPROXIMATION 

First, consider modes in the limit of the WKB approximation of 
tight- winding (i.e. local modes), where \kR\ 3> m. We caution that 
this limit does not, in fact, hold for most of the global modes seen 
in simulations, but it is nevertheless instructive. 



3.1 Instability Criteria 

The case of slow modes where the non-Keplerian part of the poten- 
tial is everywhere small is presented in detail in Tremaine (200 1| >. 
We briefly note some conclusions. Take the potential to be that of a 
central BH plus a non-Keplerian disk: 



* = $bh + $ c = • 



GAf BH 



(17) 



where /^bh ~ Mi/ 'Mbh «C 1. Further, consider "slow" modes, 
where the pattern speed fi p <C SI. We can expand the equations of 
motion in parameters of O(e) where e = Mi /Mbh <C 1. This gives 
the WKB dispersion relation (to leading order in \kR\~ 1 ) of quasi- 
Keplerian slow modes, 



uj = -cu + n GT,d \ k\ Q.~ 1 — c 2 k 2 Q.~ 1 
for a gas disk, or 

lu = vj-\-nGY>d \k\ QT T 

^m + 7rGE rf |yt|0 _1 exp(-/3\kR\) 
for a stellar disk, where we define 

q. 2 ~k 2 



(18) 



(19) 



2fi 
1 

2Q. 



2 d d~ . , 
r dr dr 2 , 



(20) 



In the dispersion relation for a stellar disk, J~ is the standard reduc- 
tion factor (Binney & Tremaine 1987), and the latter equality is a 
convenient approximation for softened gravity, with /3 ~ a z /V c ~ 
h/R (the stellar disk scale height). 

It is obvious that all terms on the right-hand side of Equa- 
tions |18|19| are real; therefore, quasi-Keplerian, "slow" m = 1 
modes are stable at this order (for a more rigorous derivation, see 
|Tremaine|200T) . 

We have, however, made a major assumption, that the disk 
mass is everywhere much less than the BH mass (and the mode is 
everywhere slow). In generality, and to second-order in the 
WKB dispersion relation can be written (for a gas disk) 

2 

(w-mfi) 2 = K 2 +(V+^)cj(l+x) (21) 

/ o m 2 \ V2 
-2nGZ d (k 2 + - r ) (1+x) (22) 
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where 



\ + {kr/m) 



l-s 

T+~s 



dlnVc 
dlnR 



(23) 



(Lau&Bertin 1978 1 



Whenever the right-hand side of Equation [22] is negative, the 
modes are unstable and grow exponentially. Take the case of inter- 
est, a global m= \ mode in a relatively cold disk. For convenience 
take the limit k = and c s = 0; Equation|22|becomes 



In- 1 ' = 2 d+*)- 2 



3-s 
T+s 



fa 



where 



ttGE M d (<R) 



tt 2 R M cnc (<R) 



(24) 



(25) 



is roughly the disk mass fraction inside R. The RHS is negative 
for fd > ( 1 + s)~ (3 — s) _ 1 . If the potential is near-Keplerian, then 
s ~ — 1 /2, so this just becomes fy > 1/10. 

In greater detail, consider the special case of a cold Mestel 
(ij — 1) disk around a BH, with an m = 1 mode and mass ratio 
enclosed in some radius y = Mj(< R) /Mbh; the full dispersion re- 
lation from Equation[22]is then 



g-l) 2 ^[l + 2,-,(|^ + D^(l + x)] 



(26) 



where for a local mode, \kR\ ^> 1 and \ 0. while for a global 
mode \kR\ ->■ but 1 + \ -> ( 7 + 6y)/( 1 + 2y). Global modes are 
formally unstable then for y > (-3 + \/i7)/4 w 0.281 (/ rf = 0.11), 
and local modes unstable for \kR\ > (1 + 2y)/y. The solutions for 
arbitrary power-law disks rj are tedious, but for global modes can 
be well approximated by y > 0.07 + 0.09 rj + 0.07 rj 1 + 0.04 rf, 

More generally, for the power-law disk+BH and the stellar 
dispersion relation, the minimum radius at which instability ap- 
pears (noting that the term \kR\ exp{— /3\kR\} is maximized for 
\kR\ = l/P) is given by 



M d {< R) 
M BH 



fie 



(2-T/)(l-/3ae(3-i7)) 



(27) 



For small /3 this is just y > 1.35/3 (1 — r)/2)~ , If the disk does not 
extend to these masses, then it will be everywhere locally stable. It 
is also immediately clear that if j3 > l/(ae(3 — rj)) (« 0.2 — 0.3 
for the interesting range of rj), then the disk is everywhere locally 
stable independent of Md/MuH (this is just Q > 1). 

This instability criteria agrees well with what is seen in sim- 
ulations; for the simulations discussed in § |T| the mode growth 
rate and maximum mode amplitudes are plotted as a function of 
Md /Menc at radii ~ lOpc near the BH radius of influence in Figure 6 
of Hopkins & Quataert (2010a) (see also Figure 12 of Hopkins & 
|Quataert|2009| >. Around these values of y or fd, rapid growth rates 
for the m = 1 mode appear at these radii. So at least at larger radii, 
mode growth is possible. 

What is the nature of these modes? Note that the right-hand 
side of Equation [26] is real; as such, to lowest order in the WKB 



' We follow jLau & Bertin|l978} keeping in-phase terms to second-order 
in because the modes of interest are global ones. One can think 

of this as accounting for the effective minimum wavenumber m from the 
azimuthal wave, and including the enhancement 1 + \ where x = T sin i 
(T = 91nr2/91nR and is the arm pitch angle), namely the leading-order 
term of the swing amplifier at this order in the WKB approximation. 



approximation, the unstable branch must correspond to an oversta- 
bility with the real part of oj, Re(o;) = Q. p = Q.. In other words, the 
system can develop fast modes that are globally unstable, where y 
is not very small. 

This suggests a picture in which the m — 1 modes first appear 
at large radii - some R C n t where M d /Mbh ~ 1 , i.e. where the poten- 
tial is only transitioning to Keplerian, and where it can be globally 
unstable. The pattern speed f2 p will simply reflect S}(R a i t ). But f 
the mode can propagate inwards at constant Q p , it will eventually 
be a slow mode, relative to the local Q. This is, in fact, what is seen 
in simulations (see §[5]below). 



3.2 Mode Structure and Propagation 

How does this occur? For now, we will remain in the WKB approx- 
imation and consider how such a mode (stable or unstable) might 
evolve. Given a mode, the wave packets themselves propagate with 
approximate group velocity v g = dw/dk — sign(^) (cj — G S)/(oj — 
f2) or 7rGE + 2c? kQ~ [ for slow modes. For a cold disk this 
is simply v s w (w — zu)R \kR\~' ; and since uj ~ Q.(R C n t ) and the 
mode is global, this is ~ V c (i?crit). The timescale for the mode to 
travel is just the dynamical time at this critical radius. 

If the mass profile is too shallow, and c s or a remains constant 
at small radii, then the wave will refract back at some Q barrier 
at some minimum radius (for constant /?, refraction occurs with 



rj < 1/2, the same criteria that Ostriker et al. 1 1992 1 show applies 



for modes in a pure fluid disk with a hard outer edge). In non- 
linear simulations, this typically leads to pile-up of inflows, grad- 
ually steepening the profile; the consequences of this for setting 
galaxy profile shapes is discussed in Hopkins & Quataert ( 2010b| >. 
Here, the mass profile is fixed; but provided the mass profile is suf- 
ficiently steep such that the RHS of Equation [T8] remains finite as 
r — > 0, then modes can propagate through to R = 0. 

Provided that the sound speed is finite, wave packets in a 
gaseous disk can propagate through the OLR to r — > 00, eventually 
becoming simple sound waves. This is discussed in Ada ms~et al.| 
( 1989) - because the waves can freely escape carrying the mode 
energy and angular momentum (and will reflect off small radii as 
above), infinite pure gaseous disks in nearly-Keplerian potentials 
do not support strong growing modes. Instead, for gaseous disks, 
mode growth is sensitive to the description of the disk edge, and if 
a "hard" edge is assumed, specifically requires efficient reflection 
of waves off the outer edge. As a consequence, it is difficult to de- 
termine the growth rates for a pure gaseous disk without some a 
priori knowledge of the edge structure, and the derivations therein 
cannot be generalized in a straightforward manner to disks with 
smooth (or thick) edges. However, for a stellar disk, the mode can- 
not propagate beyond the OLR where A = n" — in (SI — Q p )~ = 
(this acts as an effective outer edge). Refraction of the stellar waves 
off this boundary is important, and means that mode growth is pos- 
sible even when the disk extends to R 3> ^olr- Together, these con- 
straints set the dynamic range of the mode. 

Physically, how can the mode propagate, if it "begins" at larger 
radii? It is easy to see as the eccentric mode at outer radii ex- 
citing strong eccentric perturbations at smaller radii. For a slow 
m = 1 mode near the BH (where the near-Keplerian approxima- 
tion is good), the equations of motion for the perturbed velocity 
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\ = RQ^ + v r R- 



■ V0 (f> become, at this order, 

i /d<f>e 2$ f 

2 {uj n — zu) \ dr r 



{ 

2 Vr 



(28) 



(29) 



(30) 



(31) 



where & e is the "external" perturbing potential from the larger- 
scale mode. 

Consider an annulus Ri, down to which the mode has effi- 
ciently propagated, and a slightly interior radius Ro, which remains 
unperturbed. In the WKB limit the perturbing potential is domi- 
nated by the local structure - so just interior to Ri it is m $i (Ri) = 
2iv G E i (Ri ) | k\ ~ 1 . There is no local corrugation present at Ro or in- 
terior (by definition) to cancel this perturbation term. This applies 
as well in the global (non-WKB) limit; consider the limiting case 
of an element on a circular orbit inside an eccentric ring with mass 
Mring at radius Ri (with m — 1 mode amplitude \a\ defining the ec- 
centricity of the ring). The magnitude of the asymmetric term in 
the potential just inside Ri is ~ \a \ GMtmg/Ri ~ 7rGSi(i?i) Ri. If 
the disk is cold, then the local pattern speed is given directly by 
u) — w ~ 7rGEo \k\ fi -1 . Taking these estimates for the potential 
perturbation and pattern speed and applying them in Equation |28| 
we obtain 

v r = -^-\kR\- 1 QR 

U Vr Sl I! Pi-' I 

e ~ — = — \kR\ — r — - ~ \a\ . 
1 V c E 1 \kR\ 1 

For non- trivial mode amplitude \a\ ~ Ei/Eo, and a global mode 
\kR\ ~ 1, corresponding eccentricities and coherent m — 1 mode 
amplitudes are induced. This of course can then induce eccentricity 
at the next smaller annulus, and so on, allowing the perturbation 
to grow even at small R. By the same arguments, at much steeper 
slopes r\ > 1, the system will become more "stiff" against inwards 
propagation of eccentricity, but there is no strict cutoff/refraction 
(see also Zakamska & Tremaine 2004 who find the same for disks 
of planets and planetesimals). 

The key facet to note is that for an "external" driver of the 
perturbation <!>,,, the response is large independent of the local radio 
of disk to BH mass, and is linear in Thus, the mode only needs 
to be locally unstable somewhere in order to self-amplify (recall, 
our stability analysis is in the WKB limit and hence local). Self- 
gravity will grow the strength of the mode there, but the system 
will respond coherently at small radii, even where the system is 
nominally stable (a local mode there would not self-amplify). 

Finally, in both gas and stars, near the inner refraction radius 
(for r\<\ /2), an initially global mode must wind up to \k\ ~ £lc~ l 
or \k\ ~ Ocr z _1 , respectively. Near the OLR or r — > (if 77 > 1/2) 
in stars, then exp{— j3 \kR\} — > so there are the standard two 
branches: long \kR\ < /3 and short \kR\ > ft. The long branch so- 
lution corresponds to the extension of the g-modes described in 



Tremaine 



1 2001} , but for r\ > 1/2 and an OLR out at large radii 
where Mj/Mbh is not small, they do not have to have negative 
(retrograde) pattern speeds. The short branch corresponds to the 
p-modes described there. 



structure in exact solutions to the perturbed linear equations of mo- 
tion. 



4.1 Equations of Motion 

We must choose a specific disk model, so for convenience, adopt 
the power-law model in § [2] We can freely define Ro anywhere; 
so choose it such that (for the infinite disk) Ma(< Ro) = [a(2 — 
f])] Mbh - i-e. Eo = MQu/(2naRo). We can then still apply any 
disk cutoff with a in units of Ro determining the ratio Md /Mbh- We 
will consider all numerical quantities in units Mbh — G = Ro ~ 1 . 

Following Tremaine (2001), define the perturbed quantities 
as E(r) + Ei(r,r), v d (r) + vi(r,f) = rQ4> + ui(r,t)?+vi(r,t)4>, 
$(r) + $1 (r, f), and in standard fashion write the perturbation vari- 
ables in the formXi (r,r) = X\ (r, cj>,t) =X a (r) exp \i(4> — ut)]. The 
equations of motion, together with Poisson's equation, can be writ- 
ten 



i (fi — u>) E„ = - 



r dr r 

. d 2Q, 
(Q-u) — + — 
dr r 



1 

A 



Va = 



k} d (n — oj) 

2Sid? + r 

2 



dr' r'P(r,r')E a {r') 



P(r, r') = PdiiectC' - ! r) + Pindircct (r, r) 

nG , , , . TiGr 
= - — b 1/2 (r < /r > ) + — 



(32) 
(33) 

(34) 
(35) 
(36) 
(37) 



where the kernel P for the potential includes the direct (what would 
appear even in an inertial frame) and indirect components. The 
Laplace coefficient bi/2 is given by 



>1M 



cos 9 dO 



(l-2xcos6' + x 2 + /3 2 ) 1 /2 



(38) 



where /3 represents the gravitational softening." We choose this par- 
ticular softening both because it is numerically convenient, and be- 
cause it reduces exactly to the solution for a disk with scale height 
P = h/R (with p = E/2/i at |z| < h) when h < R. 

This formulation of the indirect potential follows [Murray &] 
Dermott ( 2000 1. The inclusion of the indirect potential is necessary 
because the m = 1 mode changes the center of mass, leading to mo- 
tion by the BH about that center of mass. Recall that our coordinate 
frame is centered on the BH, thus r is the vector distance to the BH, 
and the frame is rotating about the BH. 

For a given u, rj, and j3, then, the solution to these coupled 
equations determines the perturbation structure. Because we are in- 
terested in the behavior where ui/^l may be non-trivial, the equa- 
tions are non-linear in uj. However, we can combine the equation in 



4 EXACT SOLUTIONS 

The WKB approximation has been instructive. However, this does 
not allow us to survey the complete parameter space, and it is a 
suspect approximation for any global mode with small \kR\. There- 
fore, it is important to check our conclusions and examine the mode 



2 Note that when we soften the gravity for the perturbation, we also soften 
that of the unperturbed disk. This leads to a slight change in the disk poten- 
tial, and hence Q and k given in §|2]for a cold disk. However, the difference 
is very small (usually < 2%), and has no effect on our conclusions (com- 
pared to simply using f2(S) of a cold disk). 
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Figure 1. Spectrum of eigenvalues (normal mode frequencies) for m = 1 modes in the truncated power-law disk in Equation 1 1 4| around a BH, in units of 
Mbh = Rq = G = 1, as a function of the parameters f3 (the disk softening ~ h/R), 7] (the mass profile slope) and a (the outer disk scale length, which 
determines the disk-to-BH mass ratio Mj/Mbh)- Black diamonds show the pattern speed f2 p and growth rate 7 for growing modes (7 > 0). Stable modes 
(7 = 0) are shown for comparison (red circles at 7 = 10 -5 just to be displayed). Blue triangles show growing modes with Q,, < (retrograde precession), 
plotting |f2 p |. For comparison, the dashed line corresponds to 7 = \O p \. There is a large spectrum of (mostly prograde) rapidly-growing modes with 7 ~ Q p , 
and growth rates that increase in colder and relatively more massive disks. Non-zero growth rates can be present even in disks that are locally stable (Q S> 1) 
and disks with small /Mbh ~ 0.1. 



E fl , u a , and v a to obtain a single equation 



0: 



- (fi — Lj) E n 

E 



+ 



rA 



+ fn-^A)-(fi-w)] 
wl 9" 



K 



(39) 



Where $' = / dr r (dP(r, r')/dr) E (r') and $" = 
f dr'r'(d 2 P(r, r')/dr 2 )T, a (r'), and v x = dlnX/dlnr. If we 
discretize E fl into some grid in logr and apply some summation 
rule to evaluate the integrals in $, then we can write this as a linear 
operator acting on the perturbed density, 



ME fl (r') = 



(40) 



The eigenvalues ui and eigenvectors E„ represent the ex- 
act homogenous solutions to this equation. In practice, we ob- 
tain solutions follow ing the method outlined in Ap pendix B of 
Adams et al. 1 1989) . If we multiply Equation 39 by A 2 , then 
we can eliminate all occurrences of ui in the denominator and 
write the resulting matrix equation as a fifth-order equation in u>: 
(Mo + uj Mi + uj 2 M 2 + oj 3 M 3 + uj 4 M 4 + oj s M 5 ) E„ (V ) = where 
the M„ are independent of cj. With the appropriate substitutions, 
this can then be turned into a single eigenvalue equation Ms X 5 f a = 
uifa, where if M is NxN elements, M 5x 5 is 5Nx5N, and % = 



(E^ 0(uM n yE a , £>(^M„) 2 E fl , 0(cjM„) 3 E a , C(a;M„) 4 E fl ) is 
constructed from combinations of the M„, E„ and uj. It is then 
straightforward to solve for all eigenvalues and eigenvectors. 3 



4.2 Results 

4.2.1 Growth Rates 

Figure[T]illustrates the spectrum of eigenvalues of growing modes, 
in power-law disks, as a function of the softening slope {q), and 
disk cutoff radius/mass (a). For now, we simply show all eigen- 
values of the above equations that fall in the plotted range, mak- 
ing no discrimination between local or global modes. We focus on 
the growing modes, but see that in all cases there is a very large 
spectrum of stable modes; there are also decaying modes, but these 
correspond to the complex conjugate pairs of the growing modes 
shown. 

The nearly continuous lines traced out by the eigenvalues cor- 
respond to solutions that cover some finite dynamic range well in- 



3 Numerically, we typically realize this on a grid of ~ 400 — 4000 ele- 
ments evenly spaced in log r from a factor of several outside of the OLR, 
where the mode amplitude vanishes, to r ~ 10 — 5 Rq. The results are numer- 
ically converged. Normalizing such that |S B | = N, where N is the number 
of grid points, and |Row(M)| = N for each row of M, we obtain solutions 
of M S a (f) = S with typical |<5| 2 /A^ < 10" 16 . 
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Figure 2. Scaling of the largest dynamic-range growing (7 > 0) nor- 
mal modes (largest five Ti = f \a\d\ogR). Left: At fixed slope and 
a, Mj, as a function of f3. Center: As a function of slope t] (vary- 
ing a to hold Mj = Mbh)- Right: As a function of Mj (and a). In 
each column, we show the number of zero-crossings n, the mode pat- 
tern speed (black have Q p > 0, blue have Q p < 0), the growth 
rate 7, the range 1Z, the maximum ellipticity = R{/R induced by the 
mode (mode amplitude normalized so that MAX(|a|) = 1), and the max- 
imum inflow rate induced (units of Mq (G/V/o/R 3 ,) -3 / 2 ; for a BH this 
is 8400/gasMQyr^ 1 (Af B H/lO 8 M ) 3 '' 2 (Ro/4pc)" 3 / 2 ; for a protostel- 
lar/circumstellardiskO^/gasMQyr- 1 (M./Mq) 3 / 2 (tfn/lOauT 3 / 2 ). 



side the cutoff a (i.e. where the disk is still scale-free), so that the 
solution can simply be shifted in r and cj. Broadly speaking, the 
growing modes range in growth rates from 7 ~ 0.01 — 1 \£l p \, and 
the pattern speeds fl p for both stable and growing modes tend to 
fall in a similar range. We will show how this relates to the charac- 
teristic scales of the modes below. 

At fixed 77 and a (or Mj/M^h), increasing the softening fi de- 
creases the mode growth rates 7/fip, and decreases the overall frac- 
tion of unstable modes, as expected. However, growing modes ex- 
ist even for large fi > 0.25. These are global modes, where there 
is a large contribution to the instability from the indirect potential, 
which is not treated in the WKB limit. Hence, our previous conclu- 
sions regarding sufficient fi for stability should be regarded as only 
pertaining to local instabilities. Adams et alT| ( [1989^ found similar 
results, in that global instability could appear even in Q 2> 1 (ev- 
erywhere locally stable) disks. 



100 — 

80 - 
60- 
40 

20 S 

L= 



□ "a 



Bi 



100 

80 
60 
40 
20 





0.001 0.01 0.1 



1 0.0 0.5 1.0 1.5 2.0 



100 




100 





10 




10 




1 


O YVVVVVVWWV 


rf 1 
— 0.1 




0.1 






0.01 




0.01 


8 



100 
10 

1 

0.1 
0.01 



,□□! 




0.001 0.01 0.1 



1 □ 



1 0.0 0.5 1.0 1.5 2.0 

□ □r 



-1 1 



0.1 

0.001 0.01 0.1 



0.1 



-tor 



10 




□ D 















1 




nrrgi □ $p- 









0.1 








1 0.0 0.5 1.0 1.5 2.0 
10 



„ 0.001 0.01 0.1 1 0.0 0.5 1.0 1.5 2.0 

— 2r 

E -0 



5 -2 
■o 

f -4 



n = o.5 

a = 3.0 



Pan 



= 0.10 
= 1.00 



c 



0.1 
10 

1 

0.1 

2 

2 
4 




°=Bo 



-1 1 



= 0.60 
= 0.10 



0.001 0.01 0.1 

Softening p 



1 0.0 0.5 1.0 1.5 2.0 
Slope t\ 



-10 12 
log[ M„ ist /M BH ] 



Figure 3. As Figure [2] but for the largest dynamic range stable (7 = 0) 
modes. The scalings are similar - the major difference is the existence of 
stable (but not unstable) modes at Mj A/bh- 



Modes with moderate or large negative (retrograde) pattern 
speeds, Q p <C — 1 are generally not supported. Some modes with 
small retrograde pattern speeds are supported, increasingly at large 
P and/or 77. However, such modes will tend to drive outflows 
(though not in every case), so in practice may not be able to self- 
consistently build up a disk to the steep surface density profiles 
(Xi > 1/2) needed to support them in the first place. 

At fixed fi and a/Ro, modes with lower 77 (shallower disk pro- 
files) have somewhat larger growth rates, but the effect is most pro- 
nounced only as 77 — > 0. This is partly an artifact of holding a/Ro 
fixed, since doing so while decreasing 77 means that Mj/Mbh is 
larger in the lower-77 disks. But there is also a (small) real effect 
because relatively less of the disk mass is at very small radii where 
the BH dominates the potential. The range of unstable modes is 
similar, though, and we see below that these modes propagate over 
a smaller dynamic range in radius. 

At fixed 77 and fi, changing a and correspondingly Mj has the 
expected effect - at larger Mj/Mbh, growth rates (and the fraction 
of unstable modes) increases. Also, the maximum speed of the un- 
stable modes increases, as we might expect from the WKB analysis 
(since the supported mode speeds scale with Mj /Mbh in the nearly- 
Keplerian limit). There are still significant unstable modes even at 
low a and though - they appear at approximately Mj ~ 0. 1 . So 
although it is true that in the limit of vanishingly small M r i/Ms.a, 
all modes are stable, in practice given the right mass profile and 
softening, a small finite Mj/Mbh can still give interesting mode 
growth. 
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Figure 4. Example eigenmodes of the m = 1 mode in a BH+power-law disk with softened gravity. Each row shows the properties of one such mode, with 
parameters at right (disk mass profile slope S <x R^ v , softening ~ h/R, pattern speed Q p = Re(i^), mode growth rate 7 = Im(o;), and disk-to-BH mass ratio 
Mlisk)- Radii are in units of So (the radius where the enclosed disk mass ~ Mbh)- Modes are shown from an outer radius outside the OLR (where they terminate) 
to ~ 10 _3 -Rolr- Left: Mode amplitude. The absolute magnitude of the density perturbation |S n |/E, and Re(E fl /S) (the latter shows the wavenumber, i.e. 
expfi f kdr]). Center-Left: Induced radial perturbation = R a /R (equal to the eccentricity, in the near-Keplerian regime). Center-Right: Conditions for orbit- 
crossing. The perturbation radii \R l /R\X= \dR l /dR\,and \d<j>i / d(j>\. The "all" line is (|i?i/fi| 2 + \dRi/dR\ 2 + \d(j>i /d^ 1 ) 1 1 1 . If any are > 1 , there are orbit 
crossings. Red line shows the criterion from Equation|44] for where gas shocks will be induced (when > 1); note this can occur even without formal orbit 
crossings. Right: Log of the gas inflow (solid) or outflow (dotted) rate \M\ driven by the mode, in units of /gasMrjH where Oo = (GMruRq 3 ) 1 / 2 . In the 
units here, 1 corresponds to 8400M Q yf-- 1 (M B h/1O 8 M ) 3 /2 (R /4pc)~ 3 /2 f or a BH, or O^OMQyr" 1 (M,/M ) 3 / 2 (Ro/lOau)" 3 / 2 for a circumstellar 
disk. The (arbitrary) mode amplitudes are normalized so MAX(|E fl |/S) = 1. The figure shows several modes with different and 7 in a fixed system. 
The largest 7 and Q f mode (top) is local; at lower Qp (middle) the modes cover a larger dynamic range, but still have moderate \kR\; the most global mode 
(bottom) is stable, and retrograde. 
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Figure 5. Figure [4] continued. Mode structure versus n, the number of 
nodes, in the growing modes of an otherwise fixed system. From top to 
bottom, n increases from 2 to > 100. The moderate-?! modes have simi- 
lar ttp and larger 7. Despite large n, some of these modes are clearly very 
long-wavelength. 
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Figure 6. Figure [4] continued. Mode structure versus mass profile slope r\ 
at otherwise fixed properties (top-to-bottom: r) = 0.05, 0.5, 1.1, 1.85). For 
shallow slopes r\ < 1/2, the mode at fixed /3 cannot be supported at R — ¥ 0, 
and so it is localized. At larger 77, the modes carry amplitude to R — > 0. 
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Figure 7. Figure [4] continued. Mode structure versus gravitational soft- 
ening/disk scale height at otherwise fixed properties (top-to-bottom: /3 = 
0.002, 0.04, 0.1, 0.2). Modes with otherwise similar properties are more 
global in softened/thick disks, and more tightly-wound in cold disks. A 
sufficiently cold disk does not support the same global modes because the 
implied growth rate at large radii is much larger than that at small radii. 
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Figure 8. Figure [4] continued. Mode structure versus disk-to-BH 
mass ratio at otherwise fixed properties (top-to-bottom: Mi^/M^n = 
0.31, 3.2, 5.6, 316.3). At low the most unstable (local, high-i:) modes 

appear first. Higher ~ Mbh allows more global modes. Once M,^ > 
Mbh, the structure of modes that exist at small radii (the nearly-Keplerian 
regime) is similar regardless of Mji^. 
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4.2.2 General Mode Properties 

To see how the mode properties scale in detail with these choices, 
we need to select some subset of modes to study. Typically, modes 
are compared at fixed n n odcs, the number of nodes or zero-crossings. 
However, because of the large scale-free range, there are many 
modes with the same n; moreover, with the large dynamic range 
involved, the n = modes are not necessarily the "most global" 
(they might cover only a limited range in radius). Since we are in- 
terested in modes that cover a large dynamic range, we define the 
proxy 



\a(R)\d\ogR 



(41) 



where the (arbitrary) normalization of the mode amplitude is set 
such that MAX(|a(if)|) = 1 (the maximum value for which the 
solution would not imply negative densities somewhere). Roughly 
speaking, a couple times 1Z gives the number of dex over which the 
mode has a significant amplitude. 

In Figure [2] we select the five growing (7 > 0) modes with 
the largest values of 1Z, and show how their properties scale with 
/3, rj (in this case, allowing a to vary to hold Md /Mbh fixed while 
rj varies), and the disk mass Mj. In Figure [3] we show the same, 
but for the largest-7?. stable (7 = 0) modes. We show a number 
of mode properties: the number of nodes « no dcs, the pattern speed 
\^l P \, growth rate 7 (for the unstable modes), and range 1Z defined 
above. We are also interested in the effect of the perturbation on 
eccentricities and, ultimately, inflow. The magnitude of the radial 
perturbation \Ri \/R from a given mode is, for linear perturbations, 



Ri 



1 

"A 



d<S>i 
~dR 



+ 



2mO, 



R(m£l-uj) 
r a -r p )/(r a + r„) 



(42) 



Defining the eccentricity as |e| = (r a — r p )/(r a + r p ) — \R\ \/R, we 
plot the maximum jej induced by the mode over all radii, given the 
same normalization condition MAX(|a(i?) |) = 1. If the eccentric- 
ities are anywhere significant, there can be orbit crossings, hence 
shocks and dissipation. So we can calculate the induced gas inflow 
rates, using the scalings derived in Hopkins & Quataert ( 2010a), for 
the inflow rates of gas driven by shocks induced by gravitational in- 
stabilities. They show 

OTsign(57 — Q p ) 



M : 



E gas R~ Q 



V?\ 



(43) 



hdlnVc/dlnR 

where F(C, $1) = /(£)S(f, includes a weak, non-linear 

function of the induced motions from the perturbation (/(C) where 
C = dRi /BR; with 1/2 < /(C) < 2 for all C) and the phase function 
5 that depends on the relative phases of the overdensity/potential 
amplitude and radial motions/shock locations (— 1 < 5 < 1), which 
determines whether the gas is inflowing (M < 0) or outflowing 
(M > 0). We show the resulting maximum \M\, given the same nor- 
malization condition. 

We see many of the same trends as in Figure [T] in greater 
detail. With increasing ft, the growth rates (and pattern speeds, 
|fi ; ,| ~ 7]) decrease, but the dynamic range covered by the modes 
increases weakly. The latter simply comes from the fact that the 
softening suppresses local modes and increases the dynamic range 
over which the disk is in effective contact. The decrease in 7 with 
/3 is especially striking if one considers the fastest-growing modes 
(mostly local, therefore not plotted here). The maximum inflow 
rates scale weakly, decreasing at high-/3 as expected due to the in- 
creased stiffness. At high /3 we see the increased prominence of 
retrograde modes as in Figure [T] and in fact these dominate the 
high-7?. modes. Recall, at small radii the WKB dispersion relation 



reduces to Equation 19 which maximizes oj at \kR\ = 1//3, giv- 
ing for the power-law disk w mM = (e _1 /3 _1 — a(2 — rf))-K CE/QR. 
Thus, when /3 > l/(ea(2-i/)), co milK < 0, and only retrograde 
modes are supported at small radii. As a consequence, prograde 
modes, although they exist, will tend to cover a relatively smaller 
dynamic range. 

We see even at fixed Md, the trend that 7 increases with in- 
creasing slope 77. The pattern speed Q p increases with 7. The inflow 
rates increase towards steeper rj, which simply comes from the fact 
that there is proportionally more disk mass at small R. There is also 
significant change in dynamic range. From r\ w — 0.5, 1Z increases 
by a factor of 2 — 3 with 77, to a peak at r\ w 1, and then declines. 
This comes from the behavior discussed in § |3.2| wherein modes 
at 77 < 1/2 cannot propagate efficiently to R — > 0, while modes at 
77 > 1 can technically propagate but do so with declining efficiency 
(see Hopkins & Quataert 2010b). The effect is large: mode dynamic 
range goes from ~ 0.1 dex to ~ 2 dex in the unstable modes or 
0.5 dex to 5 dex in the stable modes. 

Changing the disk mass M,;/Mbh, we also obtain a number 
of interesting trends. For the chosen 77 and /3, unstable modes first 
appear at Mj/Mbh ~ 0. 1 ; the criterion depends somewhat on these 
quantities, but in general we find it agrees more or less with our 
WKB instability criterion in Equation [27] At smaller Mj/Mbh, 
there are only stable modes - these persist at arbitrarily small 
Mj /Mbh and can still drive quite significant eccentricities and in- 
flow rates, and at sufficiently small Md /Mbh their properties can be 
well-approximated by the derivations in Tremaine (2001 1. Once in- 
stability arises, there is only a weak dependence of the growth rate 
and Sip on M^/Mbh, however (it is largely determined by f2 at the 
radius where Md reaches this threshold mass). Inflow rates increase 
approximately in proportion to the disk mass. 



4.2.3 Mode Structure 



Figures |4|8| illustrate some of the eigenvectors (normal modes) for 
a range of eigenvalues uj, for representative choices of 77 and f3. 
We focus on the global modes; there are of course a wide spec- 
trum of local modes available on large scales of the self-gravitating 
disk. In each Figure, we plot the mode amplitude \a(R)\, as well 
as Re(E„/E) = \a(R)\ exp{i f kdr}, which shows the winding 
and corregations of the modes. The surface density perturbations 
generically experience sharp gradients and reflect off the OLR 
(Re(cj) = fl + k, r ~ a few), although in a some cases the E n — > 
earlier, at COR (a typical factor ~ 3 smaller radius). 

We also plot the vector radial perturbation, defined for sim- 
plicity here as just R a /R (\R a /R\ is equivalent to the standard scalar 
eccentricity in the nearly-Keplerian regime; also in this regime 
\R a /R\ ~ |v a /Vc| where the velocity perturbation vi.r = dRi/dt). 
This clearly reflects the E perturbation. 

We also show the magnitude of the radial perturbation, \Ri\/R 
(defined above) along with the magnitude of its dimensionless 
derivative, \dR\/dR\, and the same for the azimuthal perturba- 
tion 0i. If these quantities exceed unity, there are orbit crossings, 
therefore shocks and gas dissipation. Note that orbit crossings can 
occur, in principle, even if these quantities do not exceed unity - 
this is only an approximate guideline. In fact, Hopkins & Quataert 
( |2010a| > show that in a disk with finite gas sound speed c s , there can 
be gas collisions and shocks, even where there are no orbit cross- 
ings. The requirement for dissipative encounters is that gas streams 
moving in the potential of the mode be compressed by the torques 
from the collisionless component at a velocity greater than c s . This 
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condition can be written 



2-rQ 1 / 


dvi,R 


2 




1 - \dR\/dR\ k \ 


dR 


+ 


dR 



where Q — c s k/ttGY,. This depends on the gas sound speed c s , so 
is not determined in the collisionless models here, but if we assume 
the gas in these disks has c, ~ a : , then we can determine £. The 
results are plotted; at small radii especially, £ 2> 1 is common - 
for disks with finite c s ~ a z in their gas, shocks will be common 
(owing to the large ~ V c radial motions of the eccentric modes) at 
small radii, even when formal orbit crossings are rare. 

Finally we show the induced inflow/outflow rates, using the 
scaling in Equation [43] to determine the gas inflow response to the 
stellar+BH potential. Because of our choice of units, these inflow 
rates are in units of f m Mmi (GMbh/-Ko) 1/2 , where / gas =M gas /M, 
is the gas mass fraction in the disk. The inflow rates can be quite 
large for the systems of interest - given typical numbers, a rate of 
~ 10~ 4 in these plots corresponds to ~ 0.1 — IMq yr~' accretion 
rates onto a supermassive BH, or ~ 10 -5 Mq vr _1 onto a star. 

Now, consider specifically Figures [4|5| which illustrate modes 
with different frequencies and radial wavenumbers, for a fixed disk 
system. We see several of the behaviors discussed in §|5] the modes 
are global in dynamic range, but can have sizeable \kR\ and a 
large number of nodes. Since the slopes here are steeper than or 
equal to r\ — 1/2, the modes can exist down to r — \ 0. Modes with 
smaller pattern speed extend to larger radii (the OLR moves out) 
and can sustain their amplitude over a larger dynamic range, while 
modes with Re(cj) 3> 1 are not supported at moderate to large radii. 
Modes with large negative pattern speed are not seen; but those with 
moderate pattern speed are present are are among the most long- 
wavelength modes. This behavior is similar to that of the g-modes 
in |Trem aine (2001 1 - however, that work concluded that when the 
disk mass is negligible compared to the BH mass, such modes can- 
not be supported in an isolated BH+disk system; here, we find that 
moderate disk-to-BH mass ratios allow for their existence. 

The induced velocity perturbation and eccentricities are large, 
and reflect the mode density structure. Where present, orbit cross- 
ings tend to occur between the COR and OLR (R ~ Ro), and at the 
smallest radii. However, one is rarely in the strong orbit-crossing 
regime - there seems to almost be a maximum \dR\/dR\ ~ 1. This 
is easy to understand from the WKB analysis: recall, we showed 
that in the nearly-Keplerian, low disk mass limit, |v,-/V c | ~ |a|/|W?|. 
Since we are nearly Keplerian, it is trivial to show that correspond- 
ingly \R\/R\ « \v r /Vc\. And in the WKB limit, d/dR = ik, so 
\dR\/dR\ w \a\. Since \a\ < 1 always, we do not expect to see dra- 
matic orbit crossings. However, the criteria for shocks, from |Ho]>] 
|kins & Quataert (2010a ), is satisfied over a wide range of radii with 
£ S> 1. This occurs because the perturbation velocities are a fixed 
fraction ~ \a\/\kR\ of V c , and so diverge at small radii - thus even 
without formal orbit crossings, gas is being compressed in the mode 
at velocities 3> c s at small radii, generating shocks and dissipation. 

This, correspondingly, allows for the inflow rates calculated 
(wherever £ > 1). The inflow rates drop towards R — > 0, as 
seen in hydrodynamic simulations (see Figures 5 & 12 in Hop- 
|kins & Qua taert 2009). But given the choice of units here, they 
are still quite large at radii ~ 10 -3 — 10 -2 Ro, where for the 
appropriate choice of Q p and r\, they often asymptote to ap- 
proximately constant M(R). Recall again the units used: if the 
value in Figure [4] is m, then the implied Eddington ratio of 
the BH is M/Meu « 892m/ gas (M BH /1O 8 M )'/ 2 (flo/lOpc)" 3 / 2 , 
or 508 m /gas (M BH / 1 8 M ) " 1 /4 ( E / 1 5 M pc " 2 ) 3/4 . So the im- 



plied accretion rates at these radii are about Eddington, for a mode 
of near-maximal strength. This is important, since these are the 
radii where we expect the viscous disk to take over - inside such 
a radius, the ratio of the disk to BH mass is as low as ~ 10~ 5 , 
so Q becomes extremely large and we expect star formation (and 
thus the role of the collisionless component) to be inefficient. For a 
circum-stellar disk, these radii approach the stellar radius itself. 

In Figure [6] we consider how the modes depend on the disk 
mass profile slope rj. When the slope is shallow, the modes are con- 
fined to larger radii; at rj ~ 0.5 — 1.0 mode amplitudes propagate to 
R — > 0; and at much larger rj the relative amplitude can be damped 
at smaller R making orbit crossings more concentrated in R. 

In Figure[7] we compare the mode structure (at otherwise sim- 
ilar properties) to the gravitational softening or disk thickness. At 
fixed lo and other disk properties, the modes in colder disks are 
higher-/: (more tightly wound). This is also evident in the simu- 
lations in Hopkins & Quataert (2009) (see their Figures 2 & 4). 
This is because, for almost all the modes of interest, the mode has 
at least some contributions from the short-branch regime - they 
are analogous to the "p-modes" in|Tremaine (2001), both in that 
the pressure/softening effects are non-negligible, in that the char- 
acteristic modes are trailing (the kR > branch of the p-modes 
remains kR > after refraction), and in that they have positive 
(prograde) pattern speeds. Because these modes depend on some 
non-zero /3, as the modes themselves heat up the disk when they go 
non-linear (ultimately stabilizing it at some Q threshold), they can 
become more global. This has the effect (at fixed mode amplitude) 
of actually increasing the efficiency of the modes at driving large 
eccentricity and inflows to small radii, although the mode growth 
rates are lower. However, most of this difference depending on disk 
thickness is concentrated at small softening; once moderate disk 
thickness > 0.05 — 0.1 is reached, the effect of the modes is actu- 
ally fairly weakly dependent on the thickness. 

Figure [8] compares the mode structure at fixed rj and ft but 
varying Mdisk (and correspondingly a). At low Mh^/Mbh ~ 0.1 — 
0.3, the first unstable modes to appear are, unsurprisingly, the maxi- 
mally unstable modes with large \kR\ « 1//3. As Mdisk increases, the 
spectrum of unstable modes expands to include longer-wavelength 
modes, and by Mdisk > Mbh includes very global modes. Once 
Mdisk 3> Mbh, the structure at the radii we are interested in - i.e. 
in the quasi-Keplerian regime inside ~ Ro, is essentially indepen- 
dent of Mdisk (it is identical to the case of an infinite power-law 
disk). Of course, there will in this limit be other modes at larger 
radii that are simply standard disk modes, but we are not interested 
in these behaviors. 



5 COMPARISON WITH SIMULATIONS 

Our analysis thus far is restricted to the linear regime. To see 
whether our key conclusions are robust in the non-linear regime, 
with gas+stellar systems (albeit still stellar-dominated), in the pres- 
ence of inflow, star formation, feedback, and a non-trivial potential, 
we briefly compare to the mode structure in the simulations of self- 
consistently formed nuclear stellar disks from Hopkins & Quataert] 
(2010c) . 

Figure [9] illustrates this in a typical nuclear scale simulation. 
We plot the enclosed disk Mj(< R) and BH mass, mode pattern 
speed Q p (and circular speed Q), and mode amplitudes, as a func- 
tion of time in a system that is initially smooth (i.e. has no in — 1 
perturbation). The mode first appears at some radii ~ i? cr i t , where 
Mj/Mbh ~ 1, with Q,, ~ fi(/? C nt). At early times, the inner disk 
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Figure 9. Nuclear mode origins, illustrated in a case study of a single, typ- 
ical zoom-in simulation of nuclear scales showing the formation of the nu- 
clear m = 1 mode that ultimately dominates accretion. Top Left: Image of 
the (saturated) mode. Brightness shows gas surface density, color (blue to 
yellow) encodes the specific SFR. Top Right: Best-fit mode pattern speed 
Q p (solid lines) in the simulation at each of several different times. Dot- 
ted line shows the circular velocity f2(S). Modes are not plotted where 
they have no measurable amplitude. Bottom Right: Enclosed disk mass 
Mj(< R), at the same times. Horizontal dotted line shows the BH mass. 
Bottom Right: Mode amplitude \a\ in the stellar disk. The mode originates 
where Mj ~ Mgn, as a locally fast mode (Q p ~ Q(R)). It then propagates 
inwards, becoming a slow mode at smaller R. With time, the inflows gen- 
erated move the inner instability radius inward, and the slowdown of the 
mode and change of mass profile shape (from those inflows) allow the mode 
to propagate to r — > 0, as shown. The stellar mode can be sustained for very 
long times, and at later times drives the gas mode (because stars dominate 
the mass). 



profile is quite shallow (or even hollow), because no inflow has yet 
reached the center; we see the resulting cutoff in the range of the 
mode at small radius. 

Two things work to push this range inwards. First, the mode 
slows down at early times, seen in Figure [9] This occurs both via 
angular momentum exchange with the bulge and disk at somewhat 
larger (~ 100 pc) radii (a resonant process not included in our anal- 
ysis), and via direct carrying of some of the angular momentum 
in wavepackets in the gas after reflection off the inner radius above 
(through the OLR, apparent in the WKB treatment). This slowdown 
occurs while the system is in transition from the overstable growth 
phase to the nonlinear, quasi-steady state. It halts once the SI,, is 
significantly below f2(_R cr j t ); in these simulations the mode pattern 
speeds tend to stabilize at values ~ 1 — 5kms~' pc~'. The process 
is analogous to the well-studied process of bar slowdown in un- 
stable disks (although obviously with the bulge replacing the halo, 
which is dynamically irrelevant here), and we refer to those stud- 
ies for further details ( Weinberg 1985 ; Hernquist & Weinberg 1992 
Athanassoula 2003 ; Martinez- Valpuesta et al. 2006), although there 
can be a non-trivial contribution to slowdown from the motion of 
the BH itself (damping via scattering off the background stars), an 
effect not included in analytic treatment (compare e.g. Adams et al 
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1989, Shu et al. 1990). As Q p decreases, the barrier in Equation 
likewise decreases, allowing the mode to strengthen and propagate 
inwards. 

Second, the mode generates substantial gas inflows, which 
have three important effects. They lead to both gas mass avail- 
able for further inflow, and form stars, which (discussed below) 
can sustain a long-lived mode at each radius. They directly move 




Figure 10. Nuclear mode structure. The radial wavenumber k versus radius 
R, fitted to modes at times near the peak of inflow for each of our nuclear- 
scale simulations. The modes are clearly global (\kR\ ~ 1) at most radii. 
Some remain low-& to small R, but others wind up (although they can still 
propagate at finite k to R — > 0). Most wind up near ~ lOOpc, the effective 
OLR (but behavior here is complicated by interaction with the outer modes). 



the /rf ~ 0. 1 radius inwards, allowing the simple instability region 
to come closer to the BH. They also, correspondingly, steepen the 
surface density profile. Inflows will continuously increase the slope 
as long as inflow is "stalled" at any radius, and this is clearly re- 
flected in Figure [9] as the enclosed disk mass at R < R crit rises 
rapidly with time. At this point, because of significant eccentric- 
ity at all radii, v, ~ V c and Ri ~ Rq, and orbit crossings can occur at 
all radii where the mode is supported. Thus, strong inflows can be 
sustained in systems such as that in Figure [9] down to radii where 
the BH completely dominates the potential. 

The system quickly reaches a quasi-equilibrium state. The 
m — 1 mode has been excited down to R = 0; the mode pat- 
tern speeds stabilize and the mode amplitudes saturate. Figure [T0| 
shows the structure of the modes, specifically the relation be- 
tween R and kR, in different simulations at different times. For 
each simulation snapshot shown, we fit the density distribution 
over the radii plotted to an asymmetric m = 1 mode of the form 
Ei = |a(r)|£ (Y) exp{i(f r k(r)dr + ujt-(f>)}, where E is the 
aximuthally-averaged surface density and we allow \a(r)\ and 
\k(r)\ to vary as fifth-order polynomials in log(r) (the exact order 
makes little difference, though much higher order terms are noise- 
dominated). We then plot the instantaneous R and kR for each wave. 
Much of the behavior described above is evident: the modes are 
global \kR\ ~ 1 at most radii - especially their origin - but many 
wind up at small R < 1 pc and large R > lOOpc. We see both short 
and long branches appear in our simulations, but the short branch 
is more common - this appears to be set as a consequence of the 
"initial conditions" set by the earlier mode behavior, since in the 
earlier times while the mode moves to smaller R, it must wind up 
(given a shallow "pre-inflow" initial profile. At the largest radii, the 
real behavior is complicated by the interaction with the modes from 
the larger-scale inflows. At small radii, the behavior is anticipated 
by our arguments above. 
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6 DISCUSSION AND CONCLUSIONS 

Lopsided or eccentric disk (m = 1) modes in quasi-Keplerian po- 
tentials are of considerable interest for a wide range of topics in 
astrophysics. Here, we have discussed how these modes originate 
in collisionless disks, with particular focus on the case of mostly 
stellar disks surrounding a supermassive BH. We discuss the ori- 
gin, structural properties, propagation efficiencies, and resonance 
structure of such modes, comparing the simple analytic mode struc- 
ture obtained from the WKB approximation, exact numerical so- 
lutions for the global linear mode structure in idealized disk+BH 
systems, and a full treatment of hydrodynamic simulations with 
gas+stars+BHs, star formation, and stellar feedback. 

The solutions here clearly demonstrate that global, over-stable 
modes with linear growth rates 7 > are, in fact, normal modes of 
a BH+disk system and can propagate to ;• — > 0. In other words, 
our growing modes can be growing even at very small radii near 
the BH, where the enclosed disk mass is <C Mbh (here as low as 
1(T 5 Mbh). How can we reconcile this with the demonstration in 
Tremaine (2001) that slow modes in nearly- Keplerian disks are sta- 
ble? The key is that the modes arc global, and the disk is not every- 
where quasi-Keplerian (nor is the mode everywhere a slow mode). 
If we enforce a cutoff to the disk mass profile or mode amplitude 
such that the regime of the mode is entirely inside a radius where 
Md(< R) <C Mbh, then we indeed recover these previous results and 
find that growing modes are not permitted. But if the disk extends to 
sufficiently large radii such that somewhere, Mj(< R) > (/? /R) Mbh 
(the canonical self-gravity criterion), then at this radius instability 
is possible. And because the modes are global, the mode potential 
4> a at small radii r — > can still have non-trivial contributions from 
radii where the mode is not slow and the potential is not Keple- 
rian. As such, it is at least in principle possible to support growing 
modes at small radii. 

Physically, the correct interpretation is that of the eccentricity 
propagating inwards, as discussed in § |3.2| In Tremaine (2001 ), it 
is noted that the non-axisymmetric potential $„ does not have to be 
self-generating; it can be imposed as some external $ e . The disk 
response at small radii is linear in - thus, if it has some com- 
plex ui, so will the response. Here, the effective <!>,, is generated on 
large scales by the self-generating instability, where the disk is self- 
gravitating and the potential is only weakly Keplerian. The inner 
parts of the disk are stable according to the definition of |Tremaine| 
( |2001) >; their behavior here is their linear response to the growing 
$<, imposed from where the disk is self-gravitating. 

The m — 1 mode is special in this respect, in a quasi-Keplerian 
potential, because the system is near-resonance (SI m k), so that the 
excited eccentricity at small radii is independent of the local ratio of 
disk to BH mass. All other (higher-;n) modes may exist, but their 
propagation efficiencies and ability to induce large eccentricities 
(hence shocks, dissipation, and inflow) will be strongly suppressed 
by factors ~ (Md(< R) /Mbh) at small radii. 

Such modes appear as fast modes (£l p ~ (i? c rit)) at me 
Rait where the disk is moderately self-gravitating (Ma{< R) > 
{h/R)Maa for local modes; for global modes, the criterion de- 
pends on the exact mass profile but is approximately Mj(< R) > 
0.1 Mbh). In this sense they are analogous to any unstable disk 
mode. The stellar waves are bounded by an OLR at a radius typ- 
ically a factor ~a couple in radius beyond co-rotation, although 
gaseous waves can move through this resonance. For this reason, 
growth rates of modes in pure fluid disks with a hard "edge" de- 
pend on reflection off that edge, and as such are quite sensitive to 
the boundary conditions ( |Shu et al.|1990||Ostriker et al.|1992||Nel-| 



|son et al.|199 8); but in gas+stellar systems or systems with a more 
smooth outer disk, the modes can refract and self-amplify without 
much dependence on the outer disk properties. If this represents the 
"first generation" of inflows, the inner disk profile may be shallow 
or hollow, so the mode may not be supported down to r — ► and 
will reflect off of a boundary at a factor of several smaller radii. But 
the mass profile will subsequently steepen owing to these inflows, 
until propagation to R — > is possible. 

Meanwhile, non-linear effects (seen in simulations) may have 
slowed down the pattern speed by a factor of a few. These in- 
clude exchanges of angular momentum between the inner and outer 
disk/bulge/halo, and direct carrying of angular momentum in the 
induced gaseous waves (if present). The slowdown of the mode 
moves the OLR outwards and the inner radius inwards, and allows 
for more efficient propagation of the eccentricity and mass inflows. 
Once the critical slope is reached, and most of the gas mass has 
turned into stars, the mode reaches a quasi-equilibrium state. Since 
the stellar mode is fully bounded, and there is no buildup of mass at 
an intermediate radius, the pattern speed becomes nearly static, and 
the mass profile shape evolution slows down considerably. Large 
inflows are maintained to small R as long as gas is available, with a 
rate that scales roughly as M/Meaa ~ \a\, the mode amplitude. 

In terms of structural properties, the modes of interest are 
global, with \kR\ < a couple over most of their dynamic range 
(a factor of > 10 4 in radius. Indeed, in simulations, the modes 
seen are very global \kR\ -C 1 (essentially a pure elliptical/lopsided 
disk) over a significant dynamic range, but on sub-pc scales (~ 
0.1 — 0.5 pc), they can wind up into tighter spirals. Because they 
are low-fe over most of the dynamic range, their structure is rela- 
tively insensitive to the coldness of the disk, for realistic values, 
and in extremely cold systems other effects will tend to heat the 
disks to moderate values quickly. And growing modes can be sup- 
ported even in disks with large thickness h/R > 0.3. In fact, as the 
thickness (or Q) increases, the mode growth rates do decrease, but 
most normal modes become more global, increasing their ability to 
pump up eccentricities and drive large inflow rates. 

With sufficiently large mode amplitudes, there will be orbit 
crossings at many points around the elliptical orbit, but one is gen- 
erally in the marginal orbit-crossing regime (in which there is ~ 1 
orbit-crossing per orbit, near the axis of eccentricity). However, al- 
lowing for some gas in the disk with finite sound speed, the modes 
can easily drive shocks and dissipation at most radii where they 
have significant amplitude, even when there are no formal orbit 
crossings. The resulting inflow rates in this regime are treated in 
Hopkins & Quataert (2010a). Using the scalings therein, we esti- 
mate the inflow rates generated by these systems. For characteristic 
BH masses and radii of influence, these can be very large - corre- 
sponding to accretion at or near the BH Eddington limit, for moder- 
ate mode amplitudes > 0. 1 . For plausible mass profile shapes in the 
"equilibrium" range, these shocks, dissipation, and inflow rates can 
be sustained down to extremely small radii ~ 10~ 4 — 10 _3 i?BH, by 
which point the system has become a traditional viscous accretion 
disk. These modes therefore represent a viable and probably very 
important means of powering large accretion rates onto BHs, from 
~ 0.01 — lOpc scales. 

Although we have usually considered the case of modes 
around a supermassive BH, our conclusions regarding collisionless 
disks may also be applicable in a number of non-galactic regimes. 
For example, a circum-stellar disk of planetesimals, or a system 
with a large density of planets. Comparing our results to those 
in Zakamska & Tremaine (2004), who study the efficiency of ec- 
centric perturbation propagation in planetary systems, we find that 
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most of our conclusions still obtain despite the discrete nature of 
planetary systems, provided we take S to be the smeared-out (aver- 
age) surface density profile. The mechanics of excitation and char- 
acteristic frequencies are similar, and they show that the efficiency 
of eccentricity propagation decreases for 77 > 1, which we also find. 
If the mass in the collisionless (say, planetary) component of such 
disks is significant compared to the gas mass, and if the gaseous 
component can radiate efficiently when experiencing shocks, then 
our conclusions regarding inflow rates induced by these modes 
should be intact. Rescaling our predicted inflow rates to those ap- 
propriate for, say, a protostellar disk, they can be again quite large 

lO~ 4 M yr~' (R[M dlsk = M„]/100au)~ 3/2 . Of course, our 

speculation regarding the role of star formation shaping these pro- 
files on longer time scales should be modified appropriately (al- 
though there may be some analogies with planet formation in such 
disks). There may also be cases where proto-stellar disks experi- 
ence eccentric perturbations from a collisionless component (say, 
induced modes from a binary companion or passages of neighbor- 
ing stars/star-forming cores; see Kru mholz et al.|2007[ >, and these 
disks typically have mass ratios relative to the central object and 
power-law profiles in the mass range explicitly modeled here ( |Bur-| 
|rows et aT|"l 996; Hest er et al.|1996| l. Exploring these applications 
in greater detail will be an important subject of future work. 
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